subroutine calceaaa(ie, tmpe, tmpv,  i, j, k, irrepi, irrepj, irrepk, vrtin, popin) 
use mod_iop 
use mod_t
implicit none
integer, intent(in) :: vrtin(8), popin(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*8, intent(inout):: ie(*), tmpe(*), tmpv(*)
!
if(prec==8.and.igpu==0) then 
   call calceaaa_d(ie, tmpe, tmpv, i, j, k, irrepi, irrepj, irrepk, vrtin, popin)
elseif(prec==4.and.igpu==0) then 
   call calceaaa_s(ie, tmpe, tmpv, i, j, k, irrepi, irrepj, irrepk, vrtin, popin)
elseif(prec==8.and.igpu==1) then 
   call calceaaa_dg(ie, tmpe, tmpv, i, j, k, irrepi, irrepj, irrepk, vrtin, popin)
elseif(prec==4.and.igpu==1) then 
   call calceaaa_sg(ie, tmpe, tmpv, i, j, k, irrepi, irrepj, irrepk, vrtin, popin)
endif
return
end
subroutine calceaaa_d(ie, tmpe, tmpv, i, j, k, & 
              irrepi, irrepj, irrepk, vrtin, popin) 
use mod_symm
use mod_orbit
use mod_t
implicit none
integer, intent(in) :: vrtin(8), popin(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*8, intent(inout):: ie(*), tmpe(*), tmpv(*)
!
integer a, b, c, idxa, idxb, idxc, irrepi2
integer irrepc, irrepab, irrepabc, irrepa, irrepb
integer numab, inumab, ioffa, ioffb, ioffc, ioffe
integer numa, numb, numc, numi, numj, numk
integer idxi, idxj, idxk, idxe, irrepij
integer ioffo(8), ioffv(8), ioffout, idxout
integer, external :: irpdso, isymoffso, isymoffabc
real*8 dtmp
!iop = 1 for AAA iop = 2 for AAB iop = 3 for BBA iop = 4 for BBB
!
irrepij = dirprd(irrepi, irrepj)
irrepabc = dirprd(irrepij ,irrepk)
!
ioffo(1) = 0
do irrepi2 = 1, nirrep - 1
   ioffo(irrepi2 + 1) = ioffo(irrepi2) + popin(irrepi2)
enddo
!
ioffv(1) = ioffo(nirrep) + popin(nirrep)
do irrepi2 = 1, nirrep - 1
   ioffv(irrepi2 + 1) = ioffv(irrepi2) + vrtin(irrepi2)
enddo
!
idxi = ioffo(irrepi) + i 
idxj = ioffo(irrepj) + j 
idxk = ioffo(irrepk) + k
do irrepc = 1, nirrep
   irrepab = dirprd(irrepc, irrepabc)
   numab = irpdso(irrepab, vrtin, vrtin, 1)
   do irrepb = 1, nirrep
      irrepa = dirprd(irrepb, irrepab)
      numa = vrtin(irrepa)
      numb = vrtin(irrepb)
      numc = vrtin(irrepc)
      ioffa = ioffv(irrepa)
      ioffb = ioffv(irrepb)
      ioffc = ioffv(irrepc)
      ioffout = isymoffabc(irrepa, irrepb, irrepc, irrepabc, vrtin) - 1
      if(irrepa>irrepb.or.irrepb>irrepc.or.irrepa>irrepc) cycle
      if(irrepa==irrepb.and.irrepb==irrepc) then
!$omp parallel do
          do c = 3, numc
          do b = 2, c-1
          do a = 1, b-1
              idxout = ioffout + (c-1)*(c-2)*(c-3)/6 + (b-1)*(b-2)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$omp end parallel do
      elseif(irrepa<irrepb.and.irrepb==irrepc) then
!$omp parallel do
          do c = 2, numc
          do b = 1, c-1
          do a = 1, numa
             idxout = ioffout + ((c-1)*(c-2)/2+b-1)*numa +  a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$omp end parallel do
      elseif(irrepa==irrepb.and.irrepb<irrepc) then
!$omp parallel do
          do c = 1, numc
          do b = 2, numb
          do a = 1, b-1
             idxout = ioffout + (c-1)*numb*(numb-1)/2+(b-2)*(b-1)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$omp end parallel do
       elseif(irrepa<irrepb.and.irrepb<irrepc) then
!$omp parallel do
          do c = 1, numc
          do b = 1, numb
          do a = 1, numa
             idxout = ioffout + ((c-1)*numb+b-1)*numa + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$omp end parallel do
       endif
    enddo
enddo
return
end
subroutine calceaaa_s(ie, tmpe, tmpv, i, j, k, & 
              irrepi, irrepj, irrepk, vrtin, popin) 
use mod_symm
use mod_orbit
use mod_t
implicit none
integer, intent(in) :: vrtin(8), popin(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*4, intent(inout):: ie(*), tmpe(*), tmpv(*)
!
integer a, b, c, idxa, idxb, idxc, irrepi2
integer irrepc, irrepab, irrepabc, irrepa, irrepb
integer numab, inumab, ioffa, ioffb, ioffc, ioffe
integer numa, numb, numc, numi, numj, numk
integer idxi, idxj, idxk, idxe, irrepij
integer ioffo(8), ioffv(8), ioffout, idxout
integer, external :: irpdso, isymoffso, isymoffabc
real*8 dtmp
!iop = 1 for AAA iop = 2 for AAB iop = 3 for BBA iop = 4 for BBB
!
irrepij = dirprd(irrepi, irrepj)
irrepabc = dirprd(irrepij ,irrepk)
!
ioffo(1) = 0
do irrepi2 = 1, nirrep - 1
   ioffo(irrepi2 + 1) = ioffo(irrepi2) + popin(irrepi2)
enddo
!
ioffv(1) = ioffo(nirrep) + popin(nirrep)
do irrepi2 = 1, nirrep - 1
   ioffv(irrepi2 + 1) = ioffv(irrepi2) + vrtin(irrepi2)
enddo
!
idxi = ioffo(irrepi) + i 
idxj = ioffo(irrepj) + j 
idxk = ioffo(irrepk) + k
do irrepc = 1, nirrep
   irrepab = dirprd(irrepc, irrepabc)
   numab = irpdso(irrepab, vrtin, vrtin, 1)
   do irrepb = 1, nirrep
      irrepa = dirprd(irrepb, irrepab)
      numa = vrtin(irrepa)
      numb = vrtin(irrepb)
      numc = vrtin(irrepc)
      ioffa = ioffv(irrepa)
      ioffb = ioffv(irrepb)
      ioffc = ioffv(irrepc)
      ioffout = isymoffabc(irrepa, irrepb, irrepc, irrepabc, vrtin) - 1
      if(irrepa>irrepb.or.irrepb>irrepc.or.irrepa>irrepc) cycle
      if(irrepa==irrepb.and.irrepb==irrepc) then
!$omp parallel do
          do c = 3, numc
          do b = 2, c-1
          do a = 1, b-1
              idxout = ioffout + (c-1)*(c-2)*(c-3)/6 + (b-1)*(b-2)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$omp end parallel do
      elseif(irrepa<irrepb.and.irrepb==irrepc) then
!$omp parallel do
          do c = 2, numc
          do b = 1, c-1
          do a = 1, numa
             idxout = ioffout + ((c-1)*(c-2)/2+b-1)*numa +  a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$omp end parallel do
      elseif(irrepa==irrepb.and.irrepb<irrepc) then
!$omp parallel do
          do c = 1, numc
          do b = 2, numb
          do a = 1, b-1
             idxout = ioffout + (c-1)*numb*(numb-1)/2+(b-2)*(b-1)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$omp end parallel do
       elseif(irrepa<irrepb.and.irrepb<irrepc) then
!$omp parallel do
          do c = 1, numc
          do b = 1, numb
          do a = 1, numa
             idxout = ioffout + ((c-1)*numb+b-1)*numa + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$omp end parallel do
       endif
    enddo
enddo
return
end
subroutine calceaaa_dg(ie, tmpe, tmpv, i, j, k, & 
              irrepi, irrepj, irrepk, vrtin, popin) 
use mod_symm
use mod_orbit
use mod_t
implicit none
integer, intent(in) :: vrtin(8), popin(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*8, intent(inout),device :: ie(*), tmpe(*), tmpv(*)
!
integer a, b, c, idxa, idxb, idxc, irrepi2
integer irrepc, irrepab, irrepabc, irrepa, irrepb
integer numab, inumab, ioffa, ioffb, ioffc, ioffe
integer numa, numb, numc, numi, numj, numk
integer idxi, idxj, idxk, idxe, irrepij
integer ioffo(8), ioffv(8), ioffout, idxout
integer, external :: irpdso, isymoffso, isymoffabc
real*8 dtmp
!iop = 1 for AAA iop = 2 for AAB iop = 3 for BBA iop = 4 for BBB
!
irrepij = dirprd(irrepi, irrepj)
irrepabc = dirprd(irrepij ,irrepk)
!
ioffo(1) = 0
do irrepi2 = 1, nirrep - 1
   ioffo(irrepi2 + 1) = ioffo(irrepi2) + popin(irrepi2)
enddo
!
ioffv(1) = ioffo(nirrep) + popin(nirrep)
do irrepi2 = 1, nirrep - 1
   ioffv(irrepi2 + 1) = ioffv(irrepi2) + vrtin(irrepi2)
enddo
!
idxi = ioffo(irrepi) + i 
idxj = ioffo(irrepj) + j 
idxk = ioffo(irrepk) + k
do irrepc = 1, nirrep
   irrepab = dirprd(irrepc, irrepabc)
   numab = irpdso(irrepab, vrtin, vrtin, 1)
   do irrepb = 1, nirrep
      irrepa = dirprd(irrepb, irrepab)
      numa = vrtin(irrepa)
      numb = vrtin(irrepb)
      numc = vrtin(irrepc)
      ioffa = ioffv(irrepa)
      ioffb = ioffv(irrepb)
      ioffc = ioffv(irrepc)
      ioffout = isymoffabc(irrepa, irrepb, irrepc, irrepabc, vrtin) - 1
      if(irrepa>irrepb.or.irrepb>irrepc.or.irrepa>irrepc) cycle
      if(irrepa==irrepb.and.irrepb==irrepc) then
!$acc kernels
!$acc loop independent
          do c = 3, numc
!$acc loop independent
          do b = 2, c-1
!$acc loop independent
          do a = 1, b-1
              idxout = ioffout + (c-1)*(c-2)*(c-3)/6 + (b-1)*(b-2)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$acc end kernels
      elseif(irrepa<irrepb.and.irrepb==irrepc) then
!$acc kernels
!$acc loop independent
          do c = 2, numc
!$acc loop independent
          do b = 1, c-1
!$acc loop independent
          do a = 1, numa
             idxout = ioffout + ((c-1)*(c-2)/2+b-1)*numa +  a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$acc end kernels
      elseif(irrepa==irrepb.and.irrepb<irrepc) then
!$acc kernels
!$acc loop independent
          do c = 1, numc
!$acc loop independent
          do b = 2, numb
!$acc loop independent
          do a = 1, b-1
             idxout = ioffout + (c-1)*numb*(numb-1)/2+(b-2)*(b-1)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$acc end kernels
       elseif(irrepa<irrepb.and.irrepb<irrepc) then
!$acc kernels
!$acc loop independent
          do c = 1, numc
!$acc loop independent
          do b = 1, numb
!$acc loop independent
          do a = 1, numa
             idxout = ioffout + ((c-1)*numb+b-1)*numa + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$acc end kernels
       endif
    enddo
enddo
return
end
subroutine calceaaa_sg(ie, tmpe, tmpv, i, j, k, & 
              irrepi, irrepj, irrepk, vrtin, popin) 
use mod_symm
use mod_orbit
use mod_t
implicit none
integer, intent(in) :: vrtin(8), popin(8)
integer, intent(in) :: i, j, k, irrepi, irrepj, irrepk
real*4, intent(inout),device :: ie(*), tmpe(*), tmpv(*)
!
integer a, b, c, idxa, idxb, idxc, irrepi2
integer irrepc, irrepab, irrepabc, irrepa, irrepb
integer numab, inumab, ioffa, ioffb, ioffc, ioffe
integer numa, numb, numc, numi, numj, numk
integer idxi, idxj, idxk, idxe, irrepij
integer ioffo(8), ioffv(8), ioffout, idxout
integer, external :: irpdso, isymoffso, isymoffabc
real*8 dtmp
!iop = 1 for AAA iop = 2 for AAB iop = 3 for BBA iop = 4 for BBB
!
irrepij = dirprd(irrepi, irrepj)
irrepabc = dirprd(irrepij ,irrepk)
!
ioffo(1) = 0
do irrepi2 = 1, nirrep - 1
   ioffo(irrepi2 + 1) = ioffo(irrepi2) + popin(irrepi2)
enddo
!
ioffv(1) = ioffo(nirrep) + popin(nirrep)
do irrepi2 = 1, nirrep - 1
   ioffv(irrepi2 + 1) = ioffv(irrepi2) + vrtin(irrepi2)
enddo
!
idxi = ioffo(irrepi) + i 
idxj = ioffo(irrepj) + j 
idxk = ioffo(irrepk) + k
do irrepc = 1, nirrep
   irrepab = dirprd(irrepc, irrepabc)
   numab = irpdso(irrepab, vrtin, vrtin, 1)
   do irrepb = 1, nirrep
      irrepa = dirprd(irrepb, irrepab)
      numa = vrtin(irrepa)
      numb = vrtin(irrepb)
      numc = vrtin(irrepc)
      ioffa = ioffv(irrepa)
      ioffb = ioffv(irrepb)
      ioffc = ioffv(irrepc)
      ioffout = isymoffabc(irrepa, irrepb, irrepc, irrepabc, vrtin) - 1
      if(irrepa>irrepb.or.irrepb>irrepc.or.irrepa>irrepc) cycle
      if(irrepa==irrepb.and.irrepb==irrepc) then
!$acc kernels
!$acc loop independent
          do c = 3, numc
!$acc loop independent
          do b = 2, c-1
!$acc loop independent
          do a = 1, b-1
              idxout = ioffout + (c-1)*(c-2)*(c-3)/6 + (b-1)*(b-2)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$acc end kernels
      elseif(irrepa<irrepb.and.irrepb==irrepc) then
!$acc kernels
!$acc loop independent
          do c = 2, numc
!$acc loop independent
          do b = 1, c-1
!$acc loop independent
          do a = 1, numa
             idxout = ioffout + ((c-1)*(c-2)/2+b-1)*numa +  a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$acc end kernels
      elseif(irrepa==irrepb.and.irrepb<irrepc) then
!$acc kernels
!$acc loop independent
          do c = 1, numc
!$acc loop independent
          do b = 2, numb
!$acc loop independent
          do a = 1, b-1
             idxout = ioffout + (c-1)*numb*(numb-1)/2+(b-2)*(b-1)/2 + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$acc end kernels
       elseif(irrepa<irrepb.and.irrepb<irrepc) then
!$acc kernels
!$acc loop independent
          do c = 1, numc
!$acc loop independent
          do b = 1, numb
!$acc loop independent
          do a = 1, numa
             idxout = ioffout + ((c-1)*numb+b-1)*numa + a
              idxa = ioffa + a
              idxb = ioffb + b
              idxc = ioffc + c
              tmpe(idxout) = tmpv(idxout)/(ie(idxi) + ie(idxj) + ie(idxk) &
                          -ie(idxa) - ie(idxb) - ie(idxc))
          enddo
          enddo
          enddo
!$acc end kernels
       endif
    enddo
enddo
return
end
